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ABSTRACT 


A  computer  program  for  two-dimensional  slamming  of 
a  Mariner  section  is  discussed.  The  program  was  de¬ 
vised  to  examine  the  practical  applicability  of  a  previ¬ 
ously  described  numerical  method  and  to  gain  insight 
for  future  development  of  a  numerical  scheme  should 
additional  effort  appear  valuable.  A  discussion  of  this 
effort  and  the  implications  for  the  future  are  given. 
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I.  INTRODUCTION 


A  numerical  method  for  two-dimensional  water  entry  of  a 
general  convex  body  was  proposed  in  [l  ]  and  applied  in  [2]  to  a 
circular  cylinder.  The  effort  in  [2]  shows  some  hope  of  obtaining 
a  reasonable  result  if  sufficiently  small  mesh  size  is  employed. 
However,  the  method  required  further  investigation  because  of  the 
uncertainty  in  accuracy  of  the  numerical  scheme. 

Because  of  the  practical  importance  of  ship  slamming,  ex¬ 
periments  were  conducted  at  DTMB  for  a  Mariner  section.  A  com¬ 
puter  program  has  been  written  and  partially  debugged  to  aid  in 
further  development  of  the  above-mentioned  method  and  possible 
application  to  the  Mariner  section.  This  program,  although  not 
completely  checked,  yields  reasonably  good  values  of  the  velocity 
potential  in  comparison  with  analytical  values  for  a  circular  cylinder 
at  the  first  time  step. 

Based  on  the  present  knowledge,  discussions  are  also  given 
for  future  development  of  numerical  schemes  which  might  be  worth¬ 
while  in  those  cases  where  simple  fitting  theories  fail  to  yield  suf¬ 
ficient  information,  or  when  the  accuracy  of  experiments  is  in  doubt. 


2 


II.  FORMULATION  OF  THE  PROBLEM 

The  mathematical  formulation  for  the  entry  of  any  symmetric 
convex  body  into  a  two-dimensional  irrotational  incompressible  flow 
field  is  briefly  reviewed  in  the  following: 

Let  the  velocity  potential,  P  ,  in  the  moving  coordinates  x,  y 
fixed  on  the  body  (Fig.  1)  be  decomposed  into  two  parts 

P  m  P  '+  Pm  *  P '*  Vy  (1) 

where 

<Pm  ■*/</  (2) 

The  perturbation  velocity  components  u'and  V ,  which  are  not  neces¬ 
sarily  small  in  comparison  with  the  entry  velocity  V,  are  given  by 

Vp'  -  tu'+fv'  (31 

On  the  free  surface  the  pressure  is  a  constant,  say,  pm  ,  thus 

p  -  (on  the  free  surface).  (4) 

Since  the  free  surface  is  rising  with  velocity  V,  the  generalized 
Bernoulli's  equation  yields 

P-Pm  -~P  *  F  (v'**  Z'''V)*Tu'*+?(y-Vt^ 


(5) 


The  particles  on  the  free  surface  move  with  the  local  velocity,  thus 
if  the  free  surface  is  expressed  by  y-  Y(x,t)  -  0  or  by  x  -X(y,t)  ■  O 
the  kinematic  condition  is  [3] 


dY 

dt 


+  V 


o 


or 


where  the  velocity  components  are 


v  m  v '  V 

(6b)  will  not  be  employed  in  this  report.  Also  to  simplify  the  computer 
program,  the  free  surface  near  the  contact  point  will  not  be  constructed 
by  following  the  particles  as  in  [2]  ,  but  will  employ  only  equation  (6a). 

On  the  solid  boundary  (wetted  part),  the  relative  normal  velocity 
is  zero,  hence 


ax 


cos  (n,  x )  * 


(If  *  v)  com  (r>,y) 


*  Tangency  condition  of  the  free  surface  to  a  wedge  can  be  shown 
from  equation  (6b)  and  equation  (7),  assuming  conical  flow.  In 
general,  application  of  this  condition  may  be  contradictory  to 
these  equations  as  ^  may  not  be  zero. 
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where  n  is  the  normal  to  the  solid  boundary.  If  the  liquid  is  bounded 
by  a  finite  tank,  equation  (7)  should  be  applied  on  the  tank  as  well  as  on 
the  immerging  body.  If  an  infinite  domain  is  considered,  the  disturb¬ 
ance  vanishes  at  infinity,  thus 

—  u'  —  v*  *  O  at  very  large  distance 

from  the  body 

The  velocity  potential  for  ir rotational  incompressible  flow  is  governed 
by  the  Laplace  equation  V*$  m  o  and  likewise 


(8) 


O 


(9) 
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III.  NUMERICAL  SCHEME 

A.  Reduction  of  Differential  Equation 

To  avoid  extrapolation,  irregular  nets  are  constructed  with  the 
rising  free  surface.  In  general,  a  six-point  formula  is  needed  for  the 
Laplace  equation.  Let  the  five  points  surrounding  (x,  y)  be  (x  *  h^, 
y  *  kj.),  I*  i  *.  6  (see  Fig.  1).  Neglecting  the  remainder  in  the  Taylor's 
series,  one  finds 

v ** 7  *  r  j\-  ,**  ,  *  '(*  +  hi ,  y* *,• )-  +'(*,  v),  hi  k i ,  i  k*  I 

+  A-,  T  hi  ,  hi  k.  ,  t'Cx+hi ,  y  +  k{)  -  4'(X,  y)\~  0  (10a) 

where  the  points  should  be  selected  so  that 

A  •  |  ^  ,  hi  ,  r  h‘ ,  ^  ki  ,Tk?  \  +  O  (10b) 

In  general,  the  error  is  about  the  first  order  of  the  maximum  of  h^  and 
ki-  If  equal  arm  rectangular  nets  could  be  used,  the  local  truncation 
error  would  be  second  order. 

When  the  net  through  the  center  point  0  consists  of  two  perpen¬ 
dicular  lines  (h2»  h4»k^«  kj  ■  0),  equation  (10)  reduces  to  the  five- 
point  formulas  C4j 


*  I  *i  , 6 i  •••  |  is  a  determinant,  the  i'th  row  of  which  is  given. 
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where  *  4*  (*  ,  y  *k{  )  ,  d.'  «^Y*.y) 

In  the  square  net  domain, hi  ■  h3  ■  k2  ■  k4,  it  further  reduces  to 

A'  ♦  **  *  *  K  -  ♦  -  © 

This  will  be  employed  for  the  domain  bounded  by  the  initial  free  surface, 
the  line  of  symmetry  and  the  lines  x  *  L,  y«-L. 

B.  Difference  Form  for  the  Boundary  Conditions 

In  general,  the  first  derivatives  can  be  expressed  by  a  three- 
point  formula 


a  4>'  a  / 

ax  T^rr^ri 


a  / 

ay 


A;  ,  +  4(*,y) 


in  which  the  point  x,  y  is  on  the  boundary  while  the  neighboring  points 
are  taken  so  that  jh-,  kj|  #  0.  The  error  would  be  first  order  in  h^ 
and  kj.  Whenhj“0,  ki*k2,  equations  (12a,b)  reduce  to 


a  4> 

* 

ax- 

h e 

a<t>' 

ay 

k, 

(lib) 


(12a) 

(12b) 


(13a) 


(13b) 


I 


The  boundary  condition  equation  (7)  is  then  approximated  by 


cos(n,x)+  eos(»,y)  a  "  VeosCp^y) 


At  the  keel 


<b‘-  &  ss  —  V  k,  *  V  h 

By  symmetry  —  ^'across  the  line  of  symmetry,  equation  (lib) 


is  reduced  to 


2  ¥>'  *  +  <y4  -  *  -  O 

For  the  point  directly  under  the  keel,^-^*  -Vh,  h  is  the  size  of  the 
square  net  in  the  inner  domain.  Equation  (16)  reduces  to 


2  4‘, *  *4  “  3  4,  ~  Vh  »  O 
Equation  (10a)  can  be  represented  by 


A,  <fi,'+  A,  *  A,  4>s  *  A*  *  A,  <t>'  ~  A,  4>'  *  O 


where  A„  >  0. 


For  other  net  points  directly  under  the  body,  equation  (14)  can 
be  used  to  eliminate  the  boundary  point,  say  i  »  2.  One  finds,  making 


~ 4*,' ,  *■  f  A, -*•  -k%  j  -*“/>, 


in  equation  (14)  before  substituting  into  equation  (18), 
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a,  +;+ 


Aw  Cosfn jt) 
h,  cos  fay) 


Vk^  +AS  4y  *  A4  44  *  At  4>s  -  <4.  4» 


o 


or 


[*■  -  -(*. 


n,  cos  (*>,  y)  j 


(19) 


It  is  noted  that  for  convex  bodies,  gos  ^(n^y)  ^  ° 

The  infinite  domain  is  approximated  by  a  sufficiently  large  finite 
domain,  thus  the  boundary  conditions  at  the  large  distance  x  **  L  and  y«*-L 
are  Simply 


4  *  O  at  x  =  L  or  y «  -L  (20) 

Finally,  the  kinematic  and  dynamic  boundary  conditions  of  the 
free  surface  must  be  imposed.  Unfortunately  the  spray  at  the  contact 
point  cannot  be  traced  due  to  the  presence  of  a  singularity,  at  least  in 
the  present  scheme  [2j,  hence  as  in  the  analytic  studies  of  pile-up 
water  effect,  the  particle  at  the  new  contact  point  comes  from  a  neigh¬ 
boring  point  on  the  free  surface.  Using  an  explicit  scheme,  the  new 
free  surface  is  constructed  by  the  difference  approximation  to  equation  (6a) 


Y(x,t+ft) 


Y  (*,t)  + 


u  Y(*+h,t)-Y(x-h,t) 


(21) 
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The  truncation  error  of  equation  (21)  is  first  order  in  ft  but  second 
order  in  dX  •  h. 

The  dynamic  condition  on  the  free  surface  and  the  Bernoulli's 


equation  yields 


rr  =  ~[t  (v’*+  zy'v)  +  7  u'*+3(y-v*)\ 


Correct  to  first  order  in  ft  ,  the  potential  on  the  free  surface  is  given 


4>'(x,Y(x,t+rt),  t  +  tt)  3  4'(x,  Y(*,t)tt)  +  ft  +  ■—  6t 

2  t'(*  ,  r(*,t)(t)  *  (4  v'm-  i  W)st  -  v'u  U  ft-g 


where 


ay  «  Y(**h,t)-Y(x-h,t)  L  ,v  .  . 

•  s  x  a  - -  if  i  t)  exists 


3Y  s  Y(*+*,t)~  Y(x,t) 
a*  h 


otherwise 


4>'(*+h,t)-4>(x-h,t) 

2  h 


(24c ) 


r,  ^  4>'(x>  Y(*,  t+rt))  -  4>(x,  Y(x,t)) 

Y(x,  t+st)-Y(*,t) 


(24d) 
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C.  Graded  Nets 

To  reduce  the  number  of  points,  graded  nets  [5]  are  employed 
in  the  domain  below  the  free  surface.  Diagonal  cross  nets  are  used 
when  making  the  transition  to  a  coarser  subregion.  Equation  (lib)  is 
applicable  at  the  center  of  the  cross  net  by  taking  the  four  corner  points 
of  the  cross  as  neighboring  points. 

D.  Successive  Over -Relaxation  Procedure 

The  successive  over-relaxation  procedure  and  its  sufficient  con¬ 
ditions  for  convergence  are  given  in  [6],  [7],  ^ 8 J ,  and  elsewhere.  It 
seems  that  if  one  cannot  prove  the  largest  eigen  values  of  the  iteration 
matrix  are  less  than  one  for  convergence  (or  greater  than  one  for  diver¬ 
gence),  the  convergence  of  the  process  at  each  time  step  is  best  deter¬ 
mined  by  carrying  out  the  process  on  the  computer.  Let  the  difference 
equation  take  the  form 


C,  *  Ct  *  c,  +  c44>1  +  Cs  4>; 


c. 


-  < 


(25) 


where  the  residue,  R  -  0  for  the  exact  solution  of  the  difference  scheme. 
The  over-relaxation  procedure  is  to  adjust  K  to  £  such  that  the  residue 
R  becomes  (1-6>)R,  eu  being  the  over-relaxation  factor,  i.  e. 


■  if  [c,  V  *  <?«  *  c,  £  *  C4  <t>4  *  ct  4>'  +  *  (l  -eo)  & 


(26) 


where 


I  £  cu  *  2 


This  process  is  applied  successively  and  repeatedly  until  all  con¬ 
verge  by  the  criterion  that 


<  S 


The  appropriate  value  of  S  is  best  selected  empirically,  such  as  by 
comparing  the  results  for  S'  and  jg  . 


(27) 


E.  Time  Step  and  Space  Mesh  Sizes 

As  mentioned,  the  contact  point  is  a  singular  point,  the  particle 
at  the  contact  point  is  lost  in  an  infinitesimally  thin  spray.  Since  the 
new  effective  contact  point  comes  from  a  neighboring  point,  it  seems 
that  the  additional  effect  of  the  body  is  best  taken  into  account  by  taking 
a  time  step  so  that  one  of  the  neighboring  points  becomes  the  next  con¬ 
tact  point.  For  this  reason,  the  value  of  the  time  step  cannot  be  set 
arbitrarily.  Admitting  the  presence  of  a  singularity  at  the  contact  point, 

one  must  impose  a  relatively  stringent  space  mesh  size  so  that  the  mag- 
d  4>' 

nitude  of  error  in  ,  and  thus  the  pressure  p,  is  admissible.  Be¬ 

cause  the  space  mesh  size  is  limited  by  the  capacity  of  the  computer, 
one  must  therefore  set  an  approximate  value  of  St  or  a  lower  bound 
to  the  value  of  St  in  the  program.  This,  however,  does  not  insure  that 
the  truncation  error  in  time  would  be  sufficiently  small. 
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To  insure  a  good  approximation  in  the  present  formulation 
(except  possibly  in  the  infected  region  near  the  singular  contact  point), 
one  must  show  that  the  apparent  truncation  error  is  negligible  by  chang¬ 
ing  mesh  sizes  both  spacewise  and  timewise. 
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IV.  ANALYTIC  SOLUTION  FOR  THE  FIRST 
TIME  STEP, 


In  the  present  numerical  scheme,  the  effect  of  the  body  is  not 
felt  before  the  first  time  step.  Therefore,  the  free  surface  is  flat  and 
the  velocity  potential  ^  remains  zero  on  the  free  surface.  Therefore, 
the  analytic  solution  for  a  reflected  body  and  fluid  domain  can  be  used 
as  a  limiting  value  to  check  the  computer  program  and  its  accuracy. 
For  simplicity,  the  circular  cylinder  will  be  employed  for  checking 
purposes.  The  known  analytic  solution  for  flow  normal  to  obstacles 
bounded  by  two  circular  arcs  [9]  is  given  below: 

Let  all  lengths  be  nondimens ionalized  by  the  semi-wetted  width. 
In  the  complex  x  plane  (Fig.  2), 


*  *  x  +■  *  £ 


co  t  £  m 


j  (*“+  e-«) 
e‘*_  e~‘* 


(28) 


The  complex  potential,  W  is  given  by: 


W  —  4  +  i  m  rt  esc  (r>  t) 


From  equation  (28)  e »  thus 


2  m' 


(29) 


W 


2  ni 


Zni 


fefj*  -mr*  (■*)*  -«> 


(30) 
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On  the  circular  arc  (-x  side),  ~  ^  »n  being 

implicitly  defined  by  the  semi -thickness  of  the  section  which  is  cot  (*£)• 
Therefore, 


if/  -  O 


(streamline) 


(31a) 


-  2  ft 


&)*  *  (i)f 

It  is  clear  that  the  coordinate  system  is  fixed  on  the  body. 
At  infinity,  for  t  «  0 


(31b) 


d W  »*  c sc(rtt)  cot  (nt) 

-  «■  u  i  y  mm  - : - - - 1 _ _ 

di  ~  ecc't 


i.e.  V  -  Um  -  /  ,  Vm  m  o 


(32a) 


Therefore,  the  velocity  at  infinity  is  also  normalized  to  unity.  On  the 
free  surface,  <r%  *  OJ  *  ff  or  «  Oj  m  o  ,  therefore 


<&  -  O 

In  the  present  numerical  scheme  o  at  the  first  time  step  on  the 

free  surface,  x  m  0,  thus 


(32b) 


*  *  -  4  vo*  -  4  -V(y-<jr) 


(33) 
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rt  «  -/S  *  *  ~')a  “  _^4,k  O-X)1 

r,  «  -/x*  *  (%+!)*  -  -J(vr-  y)*+  ( '  ~  *)* 

By  coincidence  the  values  of  4'  obtained  on  a  desk  calculator 
agreed  excellently  with  the  values  of  ^  on  the  boundary  in  £ 2 ].  This 
fact  is  discovered  from  computer  results  (Fig.  3)  based  on  much  more 
stringent  convergent  criteria  in  the  present  computer  program.  The 
limiting  value  of  on  the  cylinder  at  the  first  time  step  is  governed 
by  both  equation  (31b)  and  equation  (33)  in  terms  of  entry  velocity,  V 


(33b) 


(33c) 


times  semi- wetted  width,  b. 
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V.  COMPARISON  OF  NUMERICAL  AND  ANALYTIC 
SOLUTION  AT  t1 

The  water  entry  of  a  circular  cylinder  at  time  t^  possesses  a 
simple  analytic  solution  equivalent  to  the  limiting  case  of  infinite  domain 
and  infinitesimal  mesh  sizes  in  the  present  numerical  formulation. 

This  yields  some  guidance  as  to  the  size  of  the  meshes,  size  of 
the  domain  and  the  magnitude  of  the  convergence  factor. 

The  same  example  as  in  [2]  and  £lOj  is  used.  A  domain  of 
10"  x  10"  is  taken  when  the  semi-wetted  width  is  1"  to  simulate  an 
infinite  fluid  domain.  The  present  computer  program  grades  the  net 
three  times.  The  inner  domain  in  this  example  is  2-1/2"  x  2-1/2" 
covered  by  1/8"  x  1/8"  square  meshes  as  in  [2]. 

At  first,  a  convergence  factor  of  6  ■  10"  ^  was  used  which  is 
very  close  to  the  convergence  criterion  used  for  the  desk  calculator 
results  reported  in  [2]  .  The  latter  criterion  makes  all  the  residues 
lie  between  *  .002  and  -.001,  when  the  diagonal  coefficient  is  nor- 
malized  to  4  and  the  value  of  ^  is  cut  off  at  10"'}.  Successive  over¬ 
relaxation  with  ■  1,  yields  considerably  lower  values  of  ^  on  the 
body  than  those  reported  in  [2]. 

Next,  over-relaxation  results  with  6  m  10"^  and  S  ■  10"^  and 
Co  »  1.25  are  obtained  for  the  same  domain  and  shown  in  Figure  3.  It 
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shows  convergence  with  respect  to  S  and  the  large  probable  error  if  a 
large  $  such  as  10*^  were  used.  Hence,  5  *  10*^  or  smaller  should 
be  assigned. 

It  is  interesting  to  note  that  a  large  value  of  03  (03  m  1.9  estimated 
on  the  basis  of  a  uniform  square  net)  was  also  tried  on  the  10"  x  10" 
domain  with  8  »  10'^.  It  was  hoped  that  a  faster  rate  of  convergence 
of  the  over-relaxation  procedure  [6j  could  be  achieved;  however,  the 
process  diverged.  This  is  due  to  the  fact  that,  while  all  other  sufficient 
conditions  for  convergence  are  satisfied,  the  graded  net  used  does  not 
satisfy  property  "A"  [*]•  This  may  be  shown  by  contradiction.  No  ad¬ 
ditional  values  of  w  have  been  tried.  Based  on  experience,  it  is  quite 
likely  03  ■  1.25  can  be  used  at  all  times.  If  the  process  tends  to  diverge 
after  a  few  cycles,  the  value  of  00  can  be  replaced  by  unity  for  further 
calculations. 

An  almost  completed  result  for  a  20"  x  20"  domain  with  S  «  10"^, 

*  1.25  shows  practically  the  same  value  of  4'  (to  almost  4  figures) 

at  the  keel.  The  computation  was  not  completed  to  conserve  machine 

time,  but  residues  spot-checked  were  already  less  than  3  x  10"^.  Thus, 

a  10"  x  10"  domain  should  be  adequate,  unless  higher  accuracy  is  desired. 

.  * 

Analytic  values  of  w  (Fig.  3)  from  equations  (31b)  and  (33)  are 
in  good  agreement  with  the  converged  numerical  results.  Remaining 
differences  are  mainly  due  to  the  mesh  size  of  the  inner  domain.  One 
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may  conclude  that  for  spacewise  convergence,  the  inner  domain  should 
have  a  fine  mesh  size  of  at  least  1/8  of  the  semi-wetted  width  at  the 
first  time  step.  For  better  accuracy,  smaller  mesh  size  is  needed. 

The  process,  however,  must  converge  both  spacewise  and  time- 
wise.  An  error  of  0.001  ft ^/sec  in  across  a  6t*. 0004  sec.  in  a 

medium  (water)  of  fi  =*  1.938  slug/ft^  results  in  an  absolute  error  Ep 
in  pressure 


£ 


P 


1.336  slug /ft* 


.  oot  ft /sec 
.0004  Sec 


46.  4 


The  maximum  pressure  at  initial  time  of  contact  from  compressible 
theory  [“]  for  an  entry  velocity  of  1  ft/sec  is 

pm  —  (QCV  *  1.93a  slug /ft*  *  4  794  ft /sec  x  /  ft /sec  -  B,Z9l**/ftm 

In  this  example,  V=  8.5  ft/sec,  pB  *  78,970  #/ft^.  Although  the 

pressure  drops  rapidly,  it  is  very  likely  that  the  maximum  pressure 

would  be  around  a  few  thousand  #/ft^  at  the  time  of  interest,  at  which 

time  the  incompressible  theory  might  be  applicable,  thus  an  error  in 

^  of  10'3  ft^/sec  appears  acceptable.  It  is  clear  that  larger  error 

in  exists  in  the  vicinity  of  the  contact  point. 

For  sufficiently  small  St,  error  in  pressure  would  be  propor- 

Ax 

tional  to  y .  Assuming  that  the  result  based  on  the  case  where 
St  ".0003227,  ( AX)  jnner  domain  *  1/8"  and  V  *  8.5  ft/sec  in  [2] 


(34) 
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is  reasonable)  then  it  would  be  desirable  to  make 


&x 
V  St 


/ 

Ti¬ 


ft 


0.3  ft/mmc  K  .0003137  fc 


-  3.790 


(35a) 


or 


AS  X  • 000322  7 

J- 

96 


O.  06  33 


(35b) 


One  must  be  cautioned  that  this  criterion  is  not  reliably  established, 
analytically  or  empirically. 
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VI.  APPLICATION  TO  MARINER  SECTION 

The  planning  for  the  application  of  the  present  method  to  a 
Mariner  section  will  be  presented.  In  view  of  the  substantial  funding 
which  would  be  required  and  uncertainty  in  the  reliability  of  the  numer¬ 
ical  procedure,  the  program  has  not  been  actually  carried  out. 

A.  Description  of  Mariner  Section 

The  offsets  of  the  Mariner  model  (1/20)  were  given  by  DTMB 
and  shown  in  Table  I.  The  shape  of  the  Mariner  model  and  pressure 
gauge  locations  are  given  in  Figure  4.  For  estimation  of  6  t  and  the 
size  of  the  domain  L,  the  time  of  contact  of  the  pressure  gauge  with 
the  water,  the  duration  of  interest  and  the  maximum  wetted  width  bmx 
corresponding  to  the  depth  of  submergence  ymx  are  shown  in  Table  II. 
These  values  are  quite  approximate,  since  the  effect  of  spray  and/or 
pile  up  of  water  is  not  included. 

B.  Size  of  Fluid  Domain 

The  number  of  net  points  would  be  proportional  to  the  square 
of  the  dimension  of  the  domain,  if  uniform  mesh  size  were  used.  Some 
reduction  in  this  number  is  achieved  by  graded  nets.  The  present  pro¬ 
gram  will  double  the  inner  domain  mesh  size  twice  while  grading  the 
domain  twice.  To  approximate  the  infinite  domain,  the  size  of  the 
domain  in  the  numerical  method  should  be  about  ten  times  the  outboard 
location  of  the  pressure  gauge  concerned,  or  larger.  If  good  accuracy 
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is  sought  at  times  near  the  end  of  the  pressure  record,  a  much  larger 
domain  would  be  required.  The  memory  capacity  of  the  computer 
(IBM  7090)  limits  the  number  of  net  points  to  about  20,000.  Since  the 
maximum  pressure  on  the  ship  bottom  moves  outboard  during  slamming 
and  the  maximum  value  at  the  keel  (at  instant  of  contact)  requires  a 
compressible  theory,  the  planning  takes  into  account  prediction  of 
pressure  at  the  gauge  2-1/4"  outboard,  as  well  as  the  gauge  at  the 
keel.  For  these  reasons,  the  size  of  the  domain  L  would  be  about  25". 

C.  Time  Interval 

It  seems  that  the  larger  the  number  of  steps  taken  before  the 
time  of  interest,  the  better  the  nonlinear  effect  would  be  included.  To 
reduce  the  machine  time,  three  steps  before  time  of  contact  of  the  pres¬ 
sure  gauge  (2-1/4"  outboard)  is  assigned  arbitrarily  as  a  compromise 
and  five  steps  after  time  of  contact  will  yield  four  values  of  pressure 
at  the  fourth  to  the  seventh  step.  Therefore,  a  minimum  of  eight  steps 
is  needed.  The  case  of  V  ■=  5.67  ft/sec  will  be  selected.  St  thus  should 
be  around  .00147/3  sec  »  .00049  sec.  Since  6t  cannot  be  arbitrarily 
assigned  in  the  present  scheme,  but  would  be  in  the  neighborhood  of 
the  assigned  6t,  this  value  will  be  set  as  0.00040  sec  for  the  case  in 


consideration. 
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D.  Space  Mesh  Size 

Assume  that  y  ^  for  the  inner  domain  should  be  less  than 
the  value  given  by  equation  (35a).  If  V  »  5.67  ft/sec,  t  m  .00040  sec, 
then  AX  »  3.80  x  5.67  x  .00040  x  12"  *  .1034".  Likewise,  if  V  *  5.67  ft/ sec 
t  *  .00049  sec,  then  Ax  •  3.80  x  5.67  x  .00049  x  12"  *  .1267". 

Therefore,  it  seems  that  one  should  take  Ax  -  1/8"  or  prefer¬ 
ably  1/  16"  or  both.  Since  b  £.  1.20",  ^  <  g*,  j  <  ~g  ,  which  also  seems 
to  be  a  reasonable  choice  based  on  results  for  the  circular  cylinder. 

E.  Approximate  Number  of  Net  Points 

Case  (a)  (*x)  inner  domain  •  1/8" 

The  size  of  the  inner  domain  is  taken  to  be  slightly  larger 
than  twice  the  gauge  location,  or  2  x  2.5"  *  5".  Use  graded  nets  com¬ 
posed  of  1/8"  mesh  size  outward  to  5"  then  1/4"  mesh  size  to  15",  1/2" 
to  25";  similarly  downward.  This  requires  6,400  net  points  plus  about 
800  added  net  points,  which  includes  the  points  added  in  the  irregular 
domain  up  to  8  steps. 

Case  (b)  (Ax)  inner  domain  =  1/16" 

Similarly,  graded  nets  are  composed  of  1/16"  mesh  size 
to  5",  1/8"  to  15"  and  1/4"  outward  to  25".  This  requires  19,600  net 
points  plus  about  1,440  added  net  points  totaling  21,040.  If  necessary, 
slight  adjustment  could  be  made  to  suit  the  memory  capacity  of  the 


computer  (IBM  7090). 


F. 


Cost  Estimate 


As  an  approximate  minimum  cost  estimate,  the  time  required 
for  convergence  of  the  over-relaxation  procedure  is  assumed  to  be 
directly  proportional  to  the  number  of  net  points;  use  is  then  made  of 
the  information  gained  from  the  "checkout"  results  for  a  circular  cyl¬ 
inder.  Based  on  an  IBM  7090  renting  for  $  500  an  hour,  7,200  net 
points  and  8  steps  would  require  $  600  x  8  or  $  4,  800, provided  the  com¬ 
puter  program  had  been  completely  checked  out.  The  case  of  21,040 
net  points  would  require  $  1,800  x  8  or  $  14,400.  It  would  be  prefer¬ 
able  to  have  both  cases,  as  a  comparison  would  indicate  the  magnitude 
of  truncation  error  present. 

It  is  concluded  that  the  cost  is  much  higher  than  the  funds  avail¬ 
able  and  the  application  of  the  present  method  to  the  Mariner  section 


will  not  be  carried  out. 
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VII.  COMPUTER  PROGRAM 

A.  General  Description 

A  computer  program  based  on  the  theory  previously  presented 
has  been  written  and  partially  debugged.  The  program  was  written  in 
the  FORTRAN  language  to  be  run  on  the  IBM  7090.  A  description  of 
the  program  and  accompanying  general  flow  diagram  will  be  presented 
in  this  section.  A  more  detailed  description  and  flow  diagram  is  found 
in  the  appendix. 

Essentially,  the  problem  may  be  thought  of  as  a  successive 
over -relaxation  problem.  The  general  scheme  of  the  program  is  to 
first  initialize  all  of  the  mesh  points;  then  a  calculation  is  made  to 
advance  from  the  initial  time,  tp  »  0,  to  what  will  be  referred  to  as  tj. 
At  tj  a  new  free  surface  is  calculated  which,  since  this  is  the  first  time 
interval,  will  be  a  horizontal  line.  On  the  horizontal  line  new  mesh 
points  are  added  to  the  domain  where  the  vertical  lines  of  the  original 
domain  cross  the  new  free  surface.  The  domain  at  this  time  will  in 
general  look  like  Figure  6.  The  velocity  potential  on  the  new  free 
surface  is  calculated  for  use  as  a  boundary  condition.  After  this  is 
completed,  successive  over -relaxation  calculations  are  again  carried 


out  on  the  new  extended  domain. 
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After  these  calculations  have  converged  sufficiently,  another 
time  increment  is  computed  to  advance  from  tj  to  t2-  Added  free 
surface  mesh  points  and  the  velocity  potential  on  the  new  free  surface 
are  calculated,  followed  by  another  round  of  successive  over -relaxation 
calculations.  This  general  procedure  is  carried  out  until  the  time  of 
interest  or  the  wetted  semi-width  of  interest  is  reached. 

B.  Description  of  Initial  Domain 

Initially,  the  domain  for  solution  is  rectangular.  To  improve 
the  accuracy  in  the  region  nearest  the  keel,  the  domain  is  subdivided 
into  a  very  fine  mesh  throughout  the  upper  left  part  of  the  domain.  A 
less  fine  mesh  is  used  over  an  outer  subregion  in  the  upper  left  area 
of  the  domain.  The  balance  of  the  domain  is  subdivided  with  a  rela¬ 
tively  coarse  mesh.  Let  h  represent  the  mesh  size  of  the  finest  mesh; 
then  the  intermediate  mesh  size  is  2h  and  the  coarse  mesh  size  is  4h. 

The  transition  from  the  finest  mesh  to  the  intermediate  mesh 
and  the  transition  from  the  intermediate  mesh  to  the  coarse  mesh  are 
done  in  the  same  manner  as  [5].  The  programming  problems  en¬ 
countered  in  computing  over  a  domain  subdivided  in  this  manner  were 
overcome  through  the  use  of  a  series  of  computed  GO  TO  statements 
in  FORTRAN  which  were  used  to  make  the  appropriate  transfers  for 
proper  computation  of  subscript  values.  The  values  of  4>  at  the  mesh 
points  of  the  original  domain  are  stored  in  an  array  sequentially  from 
left  to  right,  top  to  bottom. 
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Actually,  of  course,  equation  (lib)  is  not  used  in  the  form 
shown.  It  is  rewritten  as  a  successive  over-relaxation  equation  in 
the  form  (c.f.  equation  (26)  ) 

4>‘  <*+■)  _  (+'•**+  ^<n+»+  V  4>'Aen>)~  (to- I)#™ 

Computations  on  the  line  of  symmetry  (the  left-hand  boundary  of  the 
domain  shown  in  Figure  5)  are  based  on  equation  (16).  However,  the 
program  uses  equation  (36),  computing  subscripts  so  that0i'  is  used 
instead  of  03'  in  equation  (36). 

Points  on  the  lower  boundary  of  the  domain  represent  the 
boundary  conditions  at  infinity  and  consequently  the  value  of  0'  is 
identically  zero  on  this  boundary,  as  well  as  on  the  right-hand 
boundary. 

Referring  to  Figure  5,  a  few  examples  of  the  subscripts  which 
would  be  calculated  and  used  by  the  program  in  equation  (36)  are  given. 


4' 
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Initially  all  of  the  values  of  0’  at  the  points  in  the  original  domain 
(e.  g.  ,  Figure  5)  are  set  to  zero.  It  should  also  be  noted  that  the  con¬ 
figuration  of  points  shown  in  Figure  5  is  the  minimum  configuration 
which  the  program  will  handle. 

C.  Description  of  Mesh 

The  mesh  size  of  the  finest  mesh  and  the  extent  of  the  domain 
and  the  extent  of  the  subregions  are  specified  as  input  parameters, 
the  only  restriction  being  that  the  extents  mentioned  above  must  be 
multiples  of  the  mesh  size.  The  fine  mesh  and  the  intermediate  mesh 
will  automatically  be  square  subregions.  The  complete  domain  may 
be  rectangular  if  desired. 

D.  Progression  to  Next  Time 

In  computing  the  time  increments,  a  trial  and  error  procedure 
is  used.  The  procedure  consists  of  considering  each  of  the  points  on 
the  hull  whose  half  widths  are  multiples  of  h.  The  time  required  for 
each  point  to  become  wet  is  computed  using  the  equation 


Yb( *>  -  V/«o 

2  h 


At 


(38) 
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where  Yk  (x)  is  the  height  of  the  hull  point  above  the  original  free 
surface,  the  Yf(x)  are  points  on  the  free  surface,  and  u  and  v  are 
found  by  equations  (6c)  and  (6d).  This  equation  is  derived  from  equa¬ 
tion  (21).  Each  of  these  time  increments  is  compared  against  an  input 

\ 

parameter  which  specifies  the  assigned  time  increment  desired.  The 
computed  time  which  is  closest  to  the  time  increment  desired  is  se¬ 
lected  for  proceeding  to  the  next  step  of  the  computations. 

The  point  on  the  hull  corresponding  to  this  time  increment  be¬ 
comes  the  new  contact  point.  For  each  mesh  interval  to  the  right  of 
this  new  contact  point,  the  new  free  surface  is  computed  using  equa¬ 
tion  (21).  If  this  is  the  first  forward  step  in  time,  this  new  free  sur¬ 
face  will  be  a  horizontal  line.  For  any  succeeding  time  increment,  the 
free  surface  will  in  general  be  a  decreasing  function  from  the  contact 
point  to  the  right  boundary.  The  velocity  potential  on  the  new  free  sur¬ 
face,  which  will  be  used  as  new  boundary  conditions,  is  computed  ac¬ 
cording  to  equations  (23)  and  (24).  Here  again  the  first  step  yields  a 
special  case;  namely,  that  the  velocity  potential  on  the  new  free  surface 
is  identically  zero. 

All  mesh  points  added  by  the  processes  outlined  above  are  placed 
sequentially  in  another  array  in  the  program.  An  example  of  the  general 
numbering  scheme  is  shown  in  Figures  6  and  7. 
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it  will  be  noted  that  in  equation  (38)  it  is  necessary  to  determine 
the  vertical  distance  of  the  hull-point  under  consideration  above  the 
original  free  surface  of  the  water.  This  is  done  by  a  subroutine,  FORM, 
which  computes  the  corresponding  vertical  distances  and  also  the  cosines 
of  the  angles  which  the  normal  to  the  hull  makes  with  the  x  and  y  axes 
at  that  point.  To  use  the  main  FORTRAN  program  for  any  convex  hull 
configuration,  it  would  be  necessary  only  to  write  a  short  and  simple 
FORM  routine  for  each  shape  desired.  In  the  current  work,  subroutines 
have  been  prepared  for  both  the  circular  cylinder  and  the  Mariner  section. 

E.  Special  Equations  in  the  Initial  Domain 

Referring  to  Figure  6,  there  are  several  points  of  the  original 
mesh  which  require  special  equations.  The  point  on  the  second  row  of 
the  original  mesh  which  lies  directly  under  the  keel  (point  26)  requires 
equation  (17),  which  is  used  in  the  form 

-  -t  vh)- («,-,) <*;(n) 

(39) 

After  each  iteration,  the  value  of  the  boundary  condition  at  the  keel 
(point  1  of  the  original  domain)  is  calculated  using  equation  (15). 

Points  2  through  5  on  Figure  6,  which  are  adjacent  to  one  hull 
point,  are  computed  by  an  equation  of  the  same  general  form  as  equa¬ 
tion  (11a);  this  is  given  in  over-relaxation  form  by: 
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where  interpolation  of  the  value  of  the  velocity  potential  on  the  hull  has 
been  done  as  in  [5]  and  h*  k3»  -h3»  -k4.  This  can  also  be  obtained  by 
equations  (13a),  (13b),  (14)  and  equation  (11a)  without  interpolation. 

Points  on  the  first  row  of  the  original  domain  (points  5  through 
24  in  Figure  6)  likewise  use  an  equation  of  the  form  of  equation  (11a)  and 
given  in  over-relaxation  form  by 
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these  points. 
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F.  Relaxation  Equations  in  the  Added  Domain 

For  points  which  lie  on  the  first  row  of  the  added  domain  three 
different  equations  are  used.  The  value  of  <t>'  at  the  first  point  on  the 
first  added  row  which  is  not  a  hull  point  (e.  g.  ,  point  5  of  Figure  7)  is 
found  from  an  equation  derived  by  substituting  equation  (14)  into  equa¬ 
tion  (18)  twice  to  interpolate  out  points  2  and  3.  Written  in  over¬ 
relaxation  form,  this  is: 
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and  cosj  (n,  x)  is  the  cosine  between  the  normal  to  the  hull  and  the 
x-axis  at  the  i'th  point. 

On  the  first  added  row,  if  a  point  is  adjacent  to  only  one  hull 
point,  the  equation  is  derived  from  equation  (11a)  in  a  manner  similar 
to  the  derivation  of  equations  (40),  but  with  h*  hi  ■  -113.  Examples  of 
such  points  are  points  6  and  7  on  Figure  7  and  the  equations  are 
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At  all  other  points  on  the  first  added  row  (e.  g.  ,  points  8  through 
24  of  Figure  7),  the  value  of  <J>'  is  computed  by  the  equation 
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This  is  derived  from  equation  (11a)  using  h»  hj»  -h^  and  then  writing 
in  over -relaxation  form. 

For  points  on  the  second  and  all  succeeding  rows,  there  are  also 
three  cases  to  distinguish.  First,  points  adjacent  to  two  hull  points  will 
use  equations  (42).  An  example  of  this  case  is  point  28  on  Figure  7. 

Second,  values  of  at  points  with  one  adjacent  point  on  the  hull 
(e.  g.  ,  points  29  and  30  on  Figure  7)  are  calculated  by  an  equation  based 
on  equation  (19): 
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The  k’s  are  as  defined  in  equations  (42). 
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Third,  all  other  points  use  an  equation  based  on  equation  (10a), 
which  when  simplified  and  written  in  over-relaxation  form  is 


^  'ln*o 


where 


{*«  [*,  (*+  +  kt) +  4 
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53  *  kg  A4  kgA  k,g  (46e) 
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The  k's  are  as  defined  in  equations  (42). 

G.  Convergence 

At  each  time  step,  the  criterion  for  convergence  is  given  by 


4‘jnj\ « ff 


(47) 
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at  all  mesh  points.  At  every  tenth  iteration  all  mesh  points  are  ex¬ 
amined  to  see  if  the  value  of  the  velocity  potential  satisfies  the  con¬ 
vergence  criterion,  equation  (47).  An  attempt  has  also  been  made  to 
establish  and  include  in  the  computer  program  a  means  to  detect  di¬ 
vergence.  However,  this  method  has  been  found  to  be  applicable  only 
to  the  first  time  step.  As  has  been  pointed  out  previously,  it  is  pos¬ 
sible  to  use  a  value  of  for  the  over- relaxation  factor  which  is  so 
large  that  divergence  occurs.  It  was  planned  that  the  program  would 
compute  its  own  over -relaxation  factor  and  if  such  over -relaxation 
factor  were  found  to  cause  divergence,* calculations  would  restart 
using  an  over -relaxation  factor  of  one,  or  alternatively  a  table  of 
OJ  's  would  be  stored  and  they  would  be  used  in  decreasing  order 
should  divergence  occur. 

H.  Calculation  of  Pressure 

After  convergence  has  been  achieved  for  a  given  time  step, 
the  velocity  potential  at  all  points  on  the  hull  is  calculated  using  equa¬ 
tion  (14).  The  pressure  on  the  body  is  then  calculated  by  the  forward 
difference  form  of  equation  (5). 

»  •  -  -  ,(**"■*  -**.•»*>  -  i  (M‘.z 

vt >] 


and  printout  of  these  values  is  made. 


*  e.  g.  ,  the  maximum  residue  grows. 


(48) 
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I.  Summary 

Figure  8  is  a  simplified  flow  diagram  of  the  IBM  7090  FORTRAN 
program  which  has  been  written  and  partially  debugged.  All  parts  of 
the  program  utilized  in  making  the  calculations  for  the  first  time  step 
have  been  thoroughly  tested  and  have  produced  the  results  outlined 
elsewhere  in  this  report.  The  mechanisms  for  adding  the  second, 
third,  and  all  successive  rows  of  points  have  not  been  completely  de¬ 
bugged.  Likewise,  the  relaxation  computations  for  the  second  and  all 
succeeding  time  steps  have  not  been  thoroughly  tested. 

Nonetheless,  the  programming  approach  has  adequately  demon¬ 
strated  its  generality  and  flexibility  in  the  computations  vhich  have  been 
made  for  the  first  time  increment  for  a  variety  of  mesh  sizes,  mesh 
extents,  over-relaxation  factors,  and  hull  shapes. 


I 
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VIII.  DISCUSSION 

A.  Drawbacks  of  the  Present  Method 

The  present  numerical  formulation  predicts  the  free  surface 
shape  and  the  free  surface  boundary  values,  and  then  the  pressure  on 
the  body.  If  truncation  errors  and  round-off  errors  both  inherent  and 
direct  can  be  reduced  to  a  negligible  fraction  of  the  quantities  sought, 
the  result  would  be  more  accurate  theoretically  than  the  simple  fitting 
theory.  It  would  be  actually  more  accurate  if  the  effect  of  the  com¬ 
pressibility  and  viscosity  is  really  negligible. 

The  chief  drawback  of  the  present  numerical  scheme  (discovered 
in  [2]  )  is  the  existence  of  a  singularity  at  the  contact  point,  at  least  at 
the  first  time  step,  in  the  limit  of  infinitesimal  net  size.  Such  singu¬ 
larities  are  known  to  occur  in  incompressible  fluid  theory  when  two 
boundary  conditions  at  one  point  cannot  be  fulfilled  simultaneously. 
Nevertheless,  the  discrepancy  with  reality  in  known  cases,  such  as  the 
pressure  distribution  on  an  infinitely  thin  flat  plate  at  an  angle  of  attack, 
decreases  rapidly  away  from  the  singularity.  Thus  the  predictions  were 
valuable.  The  presence  of  a  singularity  in  the  present  scheme,  however, 
prevents  an  accurate  description  of  the  free  surface  shape,  which  is  only 
approximately  constructed  by  the  "piled  up"  water  It  is  noted  that  if  one 
elects  to  stay  away  from  approaching  the  limiting  case  of  infinitesimal 
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space  net  size,  there  would  be  no  a  priori  guidance  to  the  accuracy  of 
the  solution,  while  the  probable  error  in  pressure  in  the  other  limiting 
case  of  finite  space  mesh  size  with  infinitesimal  time  can  be  shown  to 
diverge  at  every  point  on  the  body  due  to  finite  error  in  velocity  poten¬ 
tial.  The  inability  to  describe  the  free  surface  shape  accurately  may 
lead  to  errors  of  undesirable  magnitude  in  velocity  potential  and  pres¬ 
sure.  For  this  reason,  on  one  hand,  one  cannot  definitely  expect  a  few 
per  cent  error  in  pressure,  even  with  sufficiently  small  mesh  size.  On 
the  other  hand,  there  is  a  great  risk  of  obtaining  erroneous  results  from 
insufficiently  small  space  and  time  mesh  sizes.  It  is  possible  that  the 
present  program  may  yield  an  order  of  magnitude  estimate  or  a  result 
comparable  to  flat  plate  fitting  theory,  but  the  gain,  if  any,  may  not 
justify  the  minimum  investment  required  to  complete  the  planned 
computations  for  the  Mariner  section. 

The  second  drawback  in  the  present  program  is  the  large  demand 
on  machine  capacity  and  computing  time,  which  is,  perhaps,  inherent 
with  the  present  numerical  formulation.  The  limitation  of  machine  ca¬ 
pacity  is  a  serious  obstacle  to  insuring  a  sufficiently  small  space  mesh 
size  in  the  present  program.  The  demand  of  large  computing  time  on 
the  large  computer  (IBM  7090)  requires  a  large  sum,  which  is  yet  un¬ 
justified.  It  seems  desirable  that  the  storage  requirements(for  suffi¬ 
cient  accuracy)  should  be  reduced  by  a  factor  of  about  ten  or  more  for 
further  investigation. 
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B.  Future  Development 

The  present  numerical  formulation  is  attractive  due  to  its  sim¬ 
plicity  and  is  a  reasonable  preliminary  investigation.  Future  research 
for  a  numerical  method  for  ship  slamming  should  remove  the  drawbacks 
in  the  present  procedure. 

1.  Removal  of  Contact  Point  Singularity 

In  order  to  have  a  very  accurate  limiting  solution,  the  implicit 
contact  point  singularity  should  be  removed.  If  not,  the  free  surface 
shape  may  not  be  sufficiently  accurate  since  the  spray  root  is  not  traced, 
(a)  Compressibility  effect 

The  velocity  in  a  compressible  fluid  is  always  finite,  al¬ 
though  Mach  number  can  become  infinite  for  expansion  around  a  corner. 
Therefore,  a  velocity  singularity  will  not  exist  at  the  contact  point.  Of 
course,  this  does  not  insure  accuracy  unless  the  numerical  error  is 
negligible. 

It  is  noted  that  the  probable  requirement  of  smaller  net 
size  in  the  spray,  when  it  becomes  sharp  and  curvacious,  may  be  a  great 
hindrance  to  obtaining  a  good  result  at  later  times  (times  which  may  still 
be  quite  small).  For  the  initial  period  with  supersonic  expansion  of  semi- 
wetted  width,  the  method  of  temporal  source  can  be  used,  the  free  sur¬ 
face  being  flat  and  the  boundary  condition  on  the  body  being  imposed  ap¬ 
proximately  on  the  unperturbed  free  surface  level.  During  this  period, 
an  analytic  solution  is  more  appealing  in  accuracy  and  budgeting  than  a 


numerical  method.  Beyond  this  period  the  free  surface  shape  begins  to 
pile  up  or  spray  as  disturbances  are  propagated  with  the  local  speed  of 
sound.  However,  the  duration  in  which  a  computer  program  yields  suf¬ 
ficiently  accurate  results  may  be  quite  short.  Another  factor  is  that 
the  barotropic  compressible  flow  equations  are  more  complicated  than 
those  for  incompressible  flow  due  to  the  presence  of  an  additional  vari¬ 
able,  p,  the  pressure,  which  has  to  be  solved  simultaneously  with  the 
other  variables.  Also,  the  time  step  may  have  to  be  very  small  for 
stability  reasons.  It  is  therefore  desirable  to  examine  applicability  of 
simplified  approaches  although  there  is  no  a  priori  assurance  of  an 
exceptionally  accurate  answer. 

(b)  Iterative  scheme 

In  the  conic  flow  field  for  two-dimensional  water  entry 
[ll]  ,  the  free  surface  is  tangent  to  the  body  at  all  times  (>  0),  which 
is  quite  different  from  the  initial  shape  in  the  present  method.  An  iter¬ 
ative  scheme  was  used  to  determine  the  free  surface  shape  in  [11]. 

This  suggests  that  the  actual  free  surface  shape  in  an  incompressible 
fluid  might  be  located  in  such  a  way  that  an  implicit  contact  point  singu¬ 
larity  does  not  exist.  However,  a  good  iterative  scheme  remains  to  be 
developed  and  examined.  Although  the  number  of  iterations  of  any  iter  a 
tive  scheme  for  an  accurate  free  surface  shape  is  yet  unknown,  one 
might  expect  to  spend  several  times  more  computing  time  than  that  for 
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a  similar  program  for  a  non-iterative  scheme,  as  in  the  case  if  the 
free  surface  at  each  step  is  reconstructed  based  on  an  average  velocity 
during  this  time  step. 

2.  Reduction  of  Machine  Time 

(a)  Development  of  a  very  fast  convergent  process  for 
a  system  of  algebraic  equations 

A  recent  jump  method  developed  by  Doshi*  [l2]  can 
accelerate  the  convergence  of  the  successive  over -relaxation  method 
perhaps  40%  (through  private  discussions  with  Dr.  Doshi),  although  in¬ 
vestigations  have  been  made  on  a  restricted  class  of  iteration  matrices 
with  only  real  roots  and  some  knowledge  of  the  value  of  the  negative 
roots.  However,  due  to  the  irregularity  and  the  growth  of  the  domain 
in  the  Euler  system  and  the  complexity  of  the  differential  equation  in 
the  Lagrangian  system,  an  iterative  method  much  faster  than  the  method 
employed  remains  to  be  investigated  empirically.  However,  it  cannot 
reduce  considerably  —  and  might  even  increase — the  demand  on  the 
memory  capacity  of  a  high-speed  computer. 

(b)  Reduction  of  domain  through  conformal  mapping 
For  an  infinite  net,  it  is  possible  that  the  net  can  be 

coarser  and  coarser  at  larger  distances.  There  is  a  limit  to  the  num¬ 
ber  of  gradings  of  the  domain  in  the  computer  program.  Many  more 
branching  possibilities  must  be  added  for  an  additional  grading. 

*  Senior  Research  Engineer,  Department  of  Mechanical  Sciences, 
Southwest  Research  Institute. 
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Conformal  mapping  was  successful  in  analytic  ap¬ 
proaches  to  many  incompressible  flow  problems  by  transforming  to  a 
desired  domain  in  which  simpler  solutions  can  be  found.  In  the  present 
unsteady  problem,  the  interest  in  the  gradient  of  the  potential  during  a 
small  time  interval  greatly  increased  the  need  of  smaller  mesh  size 
for  higher  accuracy.  It  may  be  worthwhile  to  investigate  the  possi¬ 
bility  of  reduction  in  net  points  after  a  conformal  mapping.  The  gain, 
if  any,  would  depend  on  physical  judgment,  the  complexity  of  the  origi¬ 
nal  domain,  and  that  of  the  governing  differential  equations. 

(c )  Lagrangian  coordinates 

The  irregular  nets  constructed  with  each  increment  of 
time  in  the  Euler  system  are  not  very  convenient  to  program.  It  is  also 
impractical  to  find  a  new  transformation  to  map  each  new  fluid  domain 
into  a  desired  domain,  unless  a  sufficiently  large  domain  is  mapped 
once  and  for  all,  in  which  the  boundary  contour  and  free  surface  at  all 
times  of  interest  are  contained.  Both  of  these  would  become  more  com¬ 
plex  in  shape,  although  the  method  has  the  advantage  that  the  Laplacian 
equation  remains  the  governing  equation.  Should  there  be  sharp  spray 
in  the  theory,  it  is  difficult  to  add  fine  local  nets  in  the  spray  at  each 
increment  of  time  in  a  computer  program. 

Therefore,  it  is  worthwhile  to  investigate  the  use  of 
Lagrangian  coordinates  which  are  fixed  with  the  initial  domain.  There 
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is  some  minor  difficulty  in  conformal  mapping  such  as  the  need  of  in¬ 
verse  transforms  or  the  calculation  of  the  position  in  the  original  domain 
corresponding  to  the  net  constructed  in  the  new  domain.  But  once  the 
derivatives  of  new  coordinates  are  found,  it  is  valid  for  all  time  and  in¬ 
dependent  of  the  ship  section.  The  over-relaxation  method  may  still  be 
convergent  although  a  rigorous  proof  might  be  difficult.  Whether  there 
is  considerable  gain  in  accuracy  and/or  computing  time  in  this  combined 
approach  remains  to  be  investigated  —  actual  computation  may  be  neces¬ 
sary  without  sufficient  knowledge  on  the  bound  of  the  derivatives  of  the 
solution  sought.  An  attempt  should  also  be  made  to  eliminate  any  implicit 
singularity,  if  at  all  possible. 


IX.  CONCLUSION 


Due  to  the  relatively  large  funding  which  would  be  required 
to  complete  the  proposed  calculations  and  due  to  the  drawbacks  of 
the  method  previously  discussed,  it  is  concluded  that  further  effort 
on  this  program  is  not  warranted.  Additional  work  to  develop  a  new 
numerical  scheme  requiring  less  machine  time  and  predicting  pres¬ 
sure  distributions  and  impulses  more  accurately  is  required.  Cer¬ 
tainly,  however,  the  insight  and  experience  gained  in  the  effort 
reported  herein  will  prove  valuable  to  those  developing  new  approaches 


and  techniques  in  the  future. 
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COMPUTER  PROGRAM 

Included  here  is  a  more  detailed  flow  diagram  for  the  computer 
program  than  the  one  given  in  Figure  8.  The  points  of  the  mesh  are 
stored  in  two  different  arrays.  One  array,  the  PHI-array,  contains 
the  points  for  the  original  domain;  the  other  array,  the  P-array  con¬ 
tains  the  points  for  the  added  domain. 

Table  III  in  Appendix  B  indicates  the  transfers  necessary  to 
compute  properly  the  subscript  values  used  in  calculating  *he  values 
at  the  mesh  points  in  the  PHI-array.  The  first  column  of  the  table  is 
the  line  number  in  the  table;  the  corresponding  program  variable  name 
is  I.  The  second  column  is  the  number  of  times  that  this  line  is  to  be 
used  in  the  calculations  before  proceeding  to  the  next  line;  the  corres¬ 
ponding  program  variable  name  is  IV.  The  next  five  columns  are  re¬ 
spectively  the  subscript  value  for  the  mesh  value  being  computed 
followed  by  the  subscript  values  for  the  four  neighbors  being  used  in 
the  calculation.  In  each  case  the  subscript  value  is  to  be  calculated 
using  the  formula  number  shown  in  the  table  entry;  the  corresponding 
program  variable  names  are  10,  II,  12,  13,  14.  The  next  column  indi¬ 
cates  the  adjustment  necessary  to  the  space  mesh  size  as  successive 
points  are  computed;  the  corresponding  program  variable  name  is  DEL. 
The  last  column  is  the  number  of  the  equation  used  for  the  calculation. 


Several  other  program  variable  names  are  used  in  the  formulas 
of  this  appendix.  These  include  G1  which  is  the  ordinate  (and  abscissa) 
of  the  extent  of  the  finest  mesh,  G2  which  is  the  ordinate  (and  abscissa) 
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of  the  extent  of  the  intermediate  mesh,  G3X  which  is  the  extent  of  the 
complete  domain  in  the  horizontal  direction,  G3Y  which  is  the  extent 
of  the  complete  domain  in  the  vertical  direction,  and  H  which  is  the 
mesh  size  for  the  finest  mesh. 

Six  integers,  NUM1,  NUM2,  NUM3X,  NUM3Y,  NUMS,  and 
NUMTX  are  defined  by  equations  A1  -A6  following,  where  XFIXF  is 


the  FORTRAN  floating  to  fixed  function: 

NUM1  =  XFIXF  (Gl/H +.  1)  (Al) 

NUM2  =  XFIXF((G2  -  Gl)/(2.  0*H)  +  .  1)  (A2) 

NUM3X  =  XFIXF  ((G3X  -  G2)/(4.  0*H)  +  .  1)  (A3) 

NUM3Y  --  XFIXF  ((G3Y  -  G2)/(4.  0*H)  +  .  1)  (A4) 

NUMS  ~  XFIXF  (G2/(2.  0*H)  +  .  1)  (A5) 

NUMTX  -  XFIXF  (G3X/(4.  0*H)  +  .  1)  (A6) 

Formulas  referred  to  in  the  table  are  given  following: 

IV  ~  NUM2-1  (A7) 

IV  -*  NUM1/2-1  (A8) 

IV  =  NUM  3X  -  1  (A9) 

IV  ~  NUM2  (A10) 

!  IV  -  NUMS -2  (All) 


j 


IV  --  NUMS/ 2-1 


(A  12) 
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IV  =  NUMTX-2 

(A  13) 

IV  =  1 

(A  14) 

IV  =  NUM1/2-2 

(A  15) 

IV  =  NUM1-1 

(A  16) 

IV  =  NUMS-1 

(A  17) 

IV  =  NUMS/2-2 

(A  18) 

IV  =  NUM2-2 

(A  19) 

IV  =  NUM3X-2 

(A20) 

IV  =  ITAB2(l)-2 

(A21) 

IV  =  NUM1  -  ITAB2(1}  +  1 

(A22) 

IOS  =  IOS  +  2 

(A23) 

10  =  1 

(A24) 

10  =  10  +  1 

(A25) 

10  =  10  +  2 

(A26) 

IOS  =  10+1 

(A27) 

11  =  10  +  1 

(A28) 

11  =  11  +  NUM1  +  NUM2  +  5 

(A29) 

IIS  =  IIS  +  2 

(A  30) 

11  =  10  -  NUM1  -  NUM2  -  NUM3X  +  1 

(A31) 

11  =  11  +  NUM3X  +  NUMS  +  4 

(A32) 

11  =  10  -  NUMS  -  NUM3X  +  1 

(A33) 

IIS  =  IOS  +  1 

(A  34) 

11  =  2 

(A  35) 
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IIS  =  IIS  +  2  (A36) 

II  =  II  +  1  (A37) 

II  =  II  +  2  (A38) 

11  =  12  +  1  (A  39) 

12  =  12  +  NUM2  +  NUM3X  +  1  (A40) 

12  =  12  -  NUM1  -  NUM2  -  NUM3X  (A41) 

12  a  12  +  NUM1  +  NUM3X  +,  1  (A42) 

12  =  12  +  NUM2  +  2  (A43) 

12  =  12  -  NUM1  -  NUM2  -  1  (A44) 

12  =  12  -  NUM2  -  NUM3X  -  2*NUM1  -  2  (A45) 

12  =  12  -  NUM1  (A46) 

I2S  =  IOS  -  2  -  NUM1  (A47) 

12  =  12  +  2  (A  48) 

I2S  =  I2S  +  2  (A  49) 

12  =  II  -  2  (A50) 

12  =  10  -  NUMS  -  2  (A51) 

12  =  10  -  NUMS  -  NUM3X  -  1  (A52) 

12  =  12  -  NUMS  (A53) 

I2S  =  IOS  -  NUMS  -  2  (A54) 

DEL  =  2.  0*DEL  (A55) 

12  =  12  -  NUM1  +  2  (A56) 

12  =  10  -  NUMI  -  NUM2  -  NUM3X  -  1  (A57) 

12  =  12  +  1  (A58) 
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12  =  12  -  NUMS  +  2 

(A59) 

12  =  10  -  NUMTX  -  1 

(A60) 

DEL  =  H 

(A61) 

12  =  10  -  2*NUMS  -  NUM3X  -  3 

(A62) 

12  =  12  -  NUMS  -  NUM1/2  -  NUM3X 

(A63) 

13  =  13  -  1 

(A64) 

13  =  14 

(A65) 

13  =  10  +  1 

(A66) 

13  =  13  -  2 

(A67) 

I3S  =  I3S  +  2 

(A68) 

13  =  10  +  NUM1/2 

(A69) 

13  =  10  +  NUMS/ 2 

(A70) 

I3S  =  IOS  -  1 

(A71) 

13  =  0 

(A72) 

13  =  13  +  1 

(A  7  3) 

13  =  13  +  2 

(A  74) 

13  =  11 

(A75) 

14  *  13  +  NUM1  +  1 

(A  76) 

14  =  14  +  NUM2  +  1 

(A77) 

14  =  14  -  NUM1  -  NUM2  -  1 

(A  78) 

14  =  14  +  NUM2  +  NUM3X 

(A  79) 

14  =  14  +  NUM1/2  +  NUM3X  +  NUM2 

(A80) 

I4S  =  IOS  +  NUM1  +  NUM2  +  NUM3X 

(A81) 
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I4S  =  I4S  +  1 

(A82) 

14  =  14  +  NUMS  +  3 

(A8  3) 

14  =  13  +  1 

(A84) 

14  •-  10  +  NUMS  +  2 

(A85) 

1.4  =  14  +  NUM3X 

(A86) 

14  =  NUMS/ 2  +  NUM3X  +  14 

P 

00 

< 

I4S  =  IOS  +  NUMS  +  NUM3X 

(A88) 

14  =  10  +  NUMTX  +  1 

(A89) 

14  =  14  +  2 

(A90) 

14  =  NUM1  +  NUM2  +  NUM3X  +  2 

(A91) 

14  =  14  +  NUM1  +  3 

(A92) 

14  =  14  +  2*NUM1  +  NUM2  +  5 

(A9  3) 

14  =  14  +  1 

(A94) 

IV  =  NUMI  -  2 

(A95) 

o 


2200 


Is  the  "first  trial"  flag  set  to  one  ? 


Is  the  time  requirement  just  com¬ 
puted  closer  to  the  time  increment 
desired  than  the  previous  computed 
time  requirement? 


Move  to  next  hull  point 


Compute  the  subscript  value  for  the 
first  point  on  the  row  being  added  and 
store  in  the  "first  point"  table 


Has  the  maximum  wetted  width  of  in¬ 
terest  been  reached? 


Add  the  time  increment  just  selected 
to  the  total  time 


Is  this  the  first  time  step? 


Is  this  the  second  time  step? 
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f 


Compute  the  velocity  potential  and 
the  height  at  each  point  of  the  added 
row  from  the  velocity  potential  at 
adjacent  points  in  the  P-array 


Calculate  the  value  of  velocity  po¬ 
tential  at  each  point  of  the  P-array 


Calculate  the  value  of  velocity  po¬ 
tential  at  each  point  of  the  PHI- 
array  (see  Table  III,  Appendix  B) 


Have  the  values  of  the  velocity  po¬ 
tential  at  all  points  converged  suf¬ 
ficiently? 


Is  this  the  first  time  step? 
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APPENDIX  B 


TABLE  I:  Offsets  of  Mariner  Model  (1/20) 


(in.) 

x  (in.) 

y  (in.) 

x  (in.) 

y  (in.) 

x  (in.) 

0 

1.20 

0 

0 

15.  00 

11.  36 

0.  10 

2.  30 

0.25 

3.  05 

16.  00 

11. 63 

0.20 

2.  87 

0.  50 

3.  95 

16.  20 

11. 69 

0.  30 

3.  31 

1.  00 

5.  02 

17.  00 

11.93 

0.40 

3.  66 

2.  00 

6.  50 

18.  00 

12.  23 

0.  50 

3.  95 

3.  00 

7.46 

19.  00 

12.  55 

0.  60 

4.21 

4.  00 

8.  10 

20.  00 

12.89 

0.70 

4.45 

5.  00 

8.  60 

21.  00 

13.  25 

0.  80 

4.68 

6.  00 

9.  01 

22.  00 

13.  62 

0.  90 

4.  87 

7.  00 

9.  34 

23.  00 

14.  00 

1.  00 

5.  02 

8.  00 

9.62 

24.  00 

14.  40 

9.  00 

9.  87 

25.  00 

14.  82 

10.  00 

10.  12 

26.  00 

15.  29 

11.  00 

10.  36 

27.  00 

15.  75 

12.  00 

10.  61 

28.  00 

16.  25 

13.  00 

10.  85 

29.  00 

16.  80 

14.  00 

11.  10 

29.70 

17.  20 

y  vertical  height 
x  horizontal  half-breadth 


Time  and  Width  of  Interest  in  Mariner 
Model  Experiments 
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TABLE  III 

Transfer  Table  for  Program 


I 

IV 

10 

11 

12 

13 

14 

DEL 

Equation 

1 

A14 

A24 

A35 

.... 

A72 

A91 

. 

2 

A21 

A25 

A37 

.... 

A73 

A94 

Abl 

40a 

3 

A22 

A25 

A37 

— 

A73 

A94 

— 

41 

4 

A14 

A25 

A37 

— 

— 

A92 

A55 

41 

5 

A14 

A25 

A  37 

— 

A74 

A94 

— 

41 

6 

A19 

A25 

A37 

— 

A73 

A94 

.... 

41 

7 

A14 

A25 

A37 

— 

— 

A93 

A55 

41 

8 

A14 

A25 

A37 

A74 

A94 

— 

41 

9 

A25 

A37 

— 

A73 

A94 

— 

41 

10 

A14 

A26 

A38 

A57 

A75 

A76 

— 

36  or  39 

11 

A14 

A25 

A37 

A58 

A64 

A94 

— 

36 

12 

A16 

A25 

A37 

A58 

A73 

A94 

— 

36 

13 

A14 

A25 

A39 

— 

A65 

A94 

— 

36 

14 

A14 

A25 

A28 

A40 

A66 

A77 

— 

36 

15 

A14 

A25 

A37 

A58 

A64 

A94 

— 

36 

16 

A95 

A25 

A37 

A58 

A73 

A94 

— 

36 

17 

A14 

A25 

A37 

A41 

— 

A92 

— 

36 

18 

A14 

A25 

A37 

A58 

A74 

A94 

— 

36 

19 

A7 

A25 

A  37 

A58 

A7  3 

A94 

— 

36 

A14 

A25 

A39 

— 

A65 

A94 

— 

36 

21 

A14 

A25 

A28 

A42 

A75 

A78 

— 

36 

22 

A14 

A25 

A37 

A58 

A64 

A94 

— 

36 

23 

Al6 

A25 

A37 

A58 

A73 

A94 

.... 

36 

24 

A14 

A25 

A39 

— 

A65 

A94 

— 

36 

25 

A14 

A25 

A28 

A43 

A75 

A79 

— 

36 

26 

A14 

A25 

A37 

A58 

A64 

A94 

.... 

36 

27 

A95 

A2  5 

A37 

A58 

A73 

A94 

— 

36 

28 

A14 

A25 

A37 

A44 

— 

A92 

— 

36 

29 

A14 

A25 

A37 

A58 

A74 

A94 

— 

36 

A19 

A25 

A37 

A58 

A73 

A94 

36 

31 

A14 

A25 

A37 

A45 

— 

A93 

— 

36 

32 

A14 

A25 

A37 

A58 

A74 

A94 

— 

36 

33 

A25 

A37 

A58 

A73 

A94 

— 

36 

34 

A14 

A25 

A29 

A46 

A7  5 

A80 

— 

36 

35 

A14 

A27, 

A26 

A34, 

A38 

A47, 

A48 

A67, 

A71 

A81, 

A94 

— 

36 

36 

A15 

A23, 

A26 

A30, 

A38 

A48, 

A49 

A68, 

A74 

A82 

A94 

— 

36 

37 

A14 

A23, 

A26 

A36, 

A37 

A48,  ' 
A49 

A68, 

A74 

A82, 

A94 

— 

36 

38 

A14 

A25 

A37 

A58 

A74 

A94 

— 

36 

39 

A19 

A25 

A37 

A58 

A73 

A94 

— 

36 
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FIGURE  I.  SIMPLIFIED  NET  CONSTRUCTION  FOR  NUMERICAL  PROCEDURE 


PERTURBATION  VELOCITY  POTENTIAL  ON  THE  BODY  <£(FT*/SEC) 


Analytical  Theory  (infinite  domain) 

Numerical  Results  (//&“ *  //<"  mesh  in  inner  Z.SH  *1.5" 
domain ) 


Symbol 

Domain 

Convergence 
Factor  S 

Relaxation 
Factor  Co 

A 

to"*  to" 

/O'* 

/ 

O 

to"*  to“ 

/O'* 

1.25 

□ 

io'*/o" 

/O'* 

1.25 

Circular  cylinder  of  S  inch  radius 
V«  9.5  ft./acc. ,  t,  •  0.0009902  sec. 


mn 


0.4  0.6 

DISTANCE  X  (IN.) 


FIGURE  3.  COMPARISON  OF  ANALYTIC  THEORY  AND  NUMERICAL  RESULTS 
AT  THE  FIRST  TIME  STEP  FOR  A  CIRCULAR  CYLINDER 


FIGURE  6.  TYPICAL  DOMAIN  WITH  FIRST  ROW  OF  POINTS  ADDED 


FIGURE  7.  TYPICAL  DOMAIN  WITH  THREE  ROWS  OF  POINTS  ADDED 
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FIGURE  8.  SIMPLIFIED  FLOW  CHART 
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